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DYNAMICS OF VARIABLE MASS SYSTEMS 


ABSTRACT 


This report presents the results of an investigation of the effects of mass loss on the 
ft ttitnHp; behavior of spinning bodies in flight. The principal goal is to determine whether 
there are circumstances under which the motion of variable mass systems can become 
unstable in the sense that their transverse angular velocities become unbounded. 

Obviously, results from a study if this kind would find immediate application in the 
aerospace field. 

The first part of this study features a complete and mathematically rigorous 
derivation of a set of equations that govern both the translational and rotational motions of 
general variable mass systems. The remainder of the study is then devoted to the 
application of the equations obtained to a systematic investigation of the effect of various 
mass loss scenarios on the dynamics of increasingly complex models of variable mass 
systems. 

It is found that mass loss can have a major impact on the dynamics of mechanical 
systems, including a possible change in the systems stability picture. Factors such as 
nozzle geometry, combustion chamber geometry, propellant’s initial shape, size and 
relative mass, and propellant location can all have important influences on the system’s 
dynamic behavior. The relative importance of these parameters on system motion are 
quantified in a way that is useful for design purposes. 
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CHAPTER 1 


INTRODUCTION 


1.1 Background 


Studies of free and forced motions of spinning rigid bodies of various geometries 
have been, and continue to be documented in great detail in the literature. Such studies 
have led to the development of important scientific instruments (gyroscopes, etc.) and to 
the concept of spin stabilization of modern spacecraft By contrast variable mass systems 
have received relatively little attention in the literature, though they play an equally 
important role in modem technology, especially in space flight. The expression “variable 
mass system,” as used in the context of this document refers to mechanical systems that 
lose and/or gain mass while in motion. Examples of such devices abound in the 
engineering literature. They include complex systems such as aircraft rockets, and moving 
robots picking up or letting go of objects, as well as simpler systems such as water 
sprinkler systems or an inflated balloon with air loss through one or more holes. 

Variable mass systems can be divided into two classes: those with continuous mass 
variation and those with discrete mass variation. Rockets, for example, fall in the 
continuous variable mass class; and, robots picking up or releasing objects, or a moving 
vehicle dropping off some of its payload in discrete chunks, fall in the discrete variable 
mass system class. Because systems with discrete mass variation can be analyzed using 
well known principles of multi-rigid-body dynamics, the focus of the analysis presented in 
this study is on systems with continuously varying mass. 



It is clear from intuition that when the net change in mass of a system, as well as its 
mass variation rate are small, it is unnecessary to account for the change in mass in the 
study of the system’s motion. For example, automobiles are in fact variable mass systems, 
yet, no one takes mass variation into account in handling and performance studies of 
ordinary automobiles. The reason is that the rate of mass variation is viewed as negligible 
— and righdy so. On the other hand, a system that undergoes a substantial change in 
mass, especially if this occurs in a short period of time, will definitely require that mass 
variation be accounted for in the study of its motion; otherwise, any predicted response of 
the system will be far removed from its true behavior. The focus of this study is on the 
determination of the impact of mass variation on the dynamic behavior of systems with 
substantial change in mass 


1.2 Early Studies of Variable Mass Systems 

Scientific study of variable mass systems has been in progress for more than two 
hundred years; so developments in this field have a long though sporadic history. The 
earliest recorded work on the dynamics of bodies with varying mass was performed in the 
18th century by Bernoulli (1738). He was then studying the forces acting on a liquid jet 
propelled hydroreactive ship — an ancient application of the principle of jet propulsion. He 
actually derived what may be referred to as the equation of motion in this special case. The 
Czech scientist and inventor, George von Buquoy (1781-1851), was the first to pose the 
general problem of the dynamics of systems with varying mass. In 1812, he obtained his 
“motion formula” for such systems, and went on to solve a large number of examples 
based on his formula, von Buquoy’s work can be said to mark the birth of the theory of 
the dynamics of systems with varying mass. In the mean time, William Moore worked out 
his mathematical theory of rocket motion in England in 1813, and, in 1819, Poisson took a 



rather modem approach and derived the equations of motion of variable mass systems 
based on Lagrange’s general formula. In their book published in 1856, Tait and Steele 
included a section on mass variation. They postulated that mass variation produced small 
continuous impacts or impulsive forces on systems, and thus resulted in continuous change 
of velocity. This work was followed, several years later, by that of Meshcherskii, whose 
work spanned the period 1897 to 1904. He essentially laid die foundation for the 
development of variable mass dynamics as a special discipline of mechanics. He devoted 
his 160 page master’s thesis to exploring a large array of issues relevant to variable mass 
dynamics — from the derivation of equations of motion to the solution of a series of 
problems in the field. All of these early investigations of variable mass systems were 
limited in one way. They were only concerned with the study of the translational motion of 
such systems. The issue of rotational motion of such systems was not addressed until the 
mid 1940’s. 

The second world war brought with it a resurgence of interest and activity in the 
dynamics of variable mass systems, mostly in connection with rocketry. At this time, 
translational motion of such systems was relatively well understood and the main focus of 
research in variable mass dynamics began to shift to the attitude motion of such systems. 
Some of the scientific giants of this new era include Rosser et al. (1947), Gantmacher and 
Levin (1947), Rankin (1949), Ellis and McArthur (1959). The equations of rotational 
motion derived by these investigators are quite similar, and have forms similar to Euler’s 
equations for rigid bodies, with extra terms that account for mass variability. Thomson 
stands out as a major contributor to this field through his book (1961) and the companion 
papers (1965, 1966). He derived several versions of the equations of motion of variable 
mass systems and his work gives a great deal of physical insight into the behavior of 
rockets. In his study of the transverse attitude motion of a non-spinning axisymmetric 
rocket, he showed that transverse rotation rate depends on the ratio-of the distance of the 
system mass center from nozzle exit, to the transverse radius of gyration of the rocket. If 



this ratio is greater than one (the most common case), the transverse angular velocity 
decreases with time; and when the ratio is less than one, the transverse angular velocity 
increases with time. Warner and Snyder (1968) brought some refinements to Thomson’s 
work and pointed out how various simplifying assumptions can lead to drastically different 
motion predictions. Meirovitch (1970) moved work on variable mass systems one step 
further by considering the impact of mass variation on variable mass rockets. 


1.3 The Star 48 Problem 

Flaws in current understanding of the dynamics of variable mass systems was 
brought to light in the early 80’ s, when several space missions with upper stages powered 
by the Star 48 solid rocket motor were observed to exhibit anomalous behavior. 
Unexpected and unexplainable rapid growth in cone angle occurred near the end of the 
motor bum. The output of a typical rate gyro mounted on one of such flights is shown in 
Fig. 1. We note from this figure that the flight is uneventful until about two thirds into the 
motor bum, when there begins an exponential growth in transverse rate, and thus in 
nutation angle. The Star 48 is the first solid rocket motor known to produce such anomaly, 
and it differs from its predecessors mainly in its much larger size and the existence of a 
submerged nozzle construction. 

The Star 48 problem sparked another flurry of investigations [ Eke (1983), Meyer 
(1983), Mingori and Yam (1986), Flandro et al. (1987), Cochran and Kang (1991)] into 
the behavior of variable mass systems, and is one of the main factors that motivated this 
work. 


In his Ph. D dissertation work, part of which was supported by this project, Wang 
(1993) [ see also Eke and Wang (1995)] modeled rocket type variable mass systems as a 
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simple cylinder of varying mass, and produced elaborate closed-form expressions that 
describe the attitude motion of such systems for various bum geometries. The study 
showed that certain propellant bum scenarios can actually cause the transverse rates of 
short and large rocket systems to diverge in a manner similar to that observed on the Star 
48 flights. 
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Fig. 1.1 Transverse Rate vs. Time for a Typical Star 48 Flight 

1.4 This Work 

This study is an extension of Wang’s work. Like Wang’s work (1993), this study 
utilizes modem mathematical tools in the study of the dynamics of variable mass systems, 
with the general intent of making valuable contributions to the field and enhancing current 
understanding of the dynamic behavior of such systems. 




Most investigators that have attempted the derivation of dynamical equations for 
variable mass systems have relied heavily on heuristics. In the next chapter, we present 
complete and mathematically rigorous derivations of both the translational and rotational 
equations of motion of variable mass systems. The remainder of the study applies the 
equations of chapter 2 to increasingly complex models of various classes of variable mass 
systems, extracting and presenting a wealth of new information on the attitude dynamics of 
such systems. 

Throughout this document, equations and figures are numbered in the form (a.b), 
where the first number a represents the chapter in which the item appears, and the second 
number b is the actual item number within the chapter. To refer to an item, the number b is 
used if the referencing is done in the same chapter in which the item appears; otherwise, the 
full label (a.b) is utilized. 



CHAPTER 2 


EQUATIONS OF MOTION 
OF VARIABLE MASS SYSTEMS 


This chapter begins with a description of the model used to characterize variable 
mass systems in this study. This model is general enough to represent a wide variety of 
physical systems that gain or lose mass while subjected to general three-dimensional 
motion. The complete equations of modem for both rotational and translational motion of 
such systems are derived using one of the methods of analytical dynamics — Kane’s 
formalism (Kane and Levinson 1985). The merit of this approach is its efficiency. It 
produces the equations of transitional and rotational motion in one mathematically rigorous 
step, and makes it possible to clarify a lot of conceptual issues in the derivation, that have 
been very difficult to do in previous work. The equations of motion that are derived are 
then compared with those obtained by Wang (1993) and others, who used the Newton- 
Euler approach. 


2.1 Model Description 

The system of interest is shown, in its most general form, in Fig. 1. It is 
determined by the closed surface, B, and its contents. The contents of B at any given 
instant can be solid (R), fluid (G), or a mixture of both. B and its contents undergo general 
three-dimensional motion in space, and matter can flow continuously in and out of B 



during this motion. For example, parts of R can “dissolve” into G by combustion or other 
processes; and, some of such products of combustion can then flow across the boundary 
B. At any given instant of time, only the surface B and whatever happens to be inside it at 
the instant constitute "the system" for that instant. Thus, the system under consideration 
here evolves continuously, both as regards its location in space, and its material 
constitution - and hence its mass. We will use the symbol S to designate this system. 



Fig. 2.1 


Variable Mass System 



The derivation of equations of motion for a system such as the one described above 
is not as straightforward as it would normally be for a system of particles and/or rigid 
bodies of constant mass. The reason is that the basic principles of dynamics, such as 
Newton-Euler equations and Lagrange’s equations are only valid when applied to a definite 
set of particles or rigid bodies. Two choices are then possible at this point One, is to seek 
or develop new formalisms that would be valid for variable mass systems; another solution 
is to model variable mass systems in a way that allows them to be viewed as constant mass 
systems, and thus make them amenable to treatment by existing principles of dynamics. 
This latter approach is the choice adopted in this study. 


2.2 Equation Formulation Strategy 

Before formulating the dynamical equations for the system of Fig. 1, we start by 
temporarily restricting the system in some important ways. First, we assume that the closed 
boundary, B, of the system maintains a constant shape, and thus encloses a region of 
constant volume at all times. In other words, B is taken to be a rigid massless shell. As 
further help in the equation derivation process, we introduce the concept of constant mass 
systems associated with the variable mass system under study. 

We consider once more the system as shown in Fig. 1, keeping in mind that the 
outer shell is now of constant shape. At some instant of time, t,, there is a definite set of 
material particles inside B. Let us assume that this set of particles is contained in a closed 
elastic container, £ lt that is identical to B at time t,. In fact, we take the viewpoint that ]£, 
has always enclosed the exact particles that ended up in B at time t,, and that B, will 
continue to delimit these particles. Obviously, subsequent to time tj, the shape of Bi will 
deviate from that of B if it is to continue to delimit the particles that were in B at time t,. 



since some (or even all) of these particles may have exited B. Similarly, prior to time t,, the 
shape of B.i was quite different from that of B, since only some, or maybe none, of the 
particles inside then were also inside B. The shape of fi, is thus seen to vary with time, 
becoming identical to, and containing the same amount of matter as B at time t,. We note 
however, that and its contents maintain the same mass at all times. We shall represent fi, 
and its contents with the symbol Sp which will be referred to henceforth as the constant 
mass system associated with the variable mass system S at time t t . Similarly, we can define 
S. 2 » S. 3 » etc., as constant mass systems associated with S at times t 2 , t 3 , etc. Furthermore, 
we assume that there exists a special subset of the particles of S that remains within B 
throughout the interval of time of interest in this study. In fact, this set is further assumed 
to constitute a rigid body, R, that is rigidly connected to B. 

The equations of motion of any one of the constant mass systems described above 
can be formulated using any of the classical approaches (Newton's Second Law, 

Lagrange's Equations, Kane's Equations, etc.) since each system is of constant mass. 
Suppose, for example, that Kane's formalism is applied to the constant mass system Sk and 
yields 


( Fr )k + [ Fr )k~ 0 (r " 1,2 ( 2 . 1 ) 

where F r is the generalized active force on the system, F* is the generalized inertia force, 
n k is the number of degrees of freedom of Sk> and the subscript k simply indicates that we 
are dealing with the system £*. Assuming that the motion of the fluid particles of Sk relative 
to the rigid part of this constant mass system is known, Sk has six degrees of freedom — 
n k = 6. 


Eqs. (1) have the explicit form 
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* * 
Q\y 



r ~1j2|>*«i 6 


( 2 . 2 ) 


where t is time, <?,*and u r k are, respectively, the r - th generalized coordinate and generalized 
speed of the constant mass system S*. In other words, the time derivative of a given 
generalized speed of £*, is a function of generalized speeds and generalized coordinates of 
&, and possibly time. Eqs. (2) can be supplemented with kinematical equations of the form 


-t k ( k k k k \ 
Qr ~ 8r I Mg, Q\, ,Q(y> tj 


r=l,2,..., 6 


(2.3) 


Eqs. (2) and (3) constitute the equations of motion of S*, and can be written as 


y k r = h k r (yi.-o{2* r ) r = 1 ’ 2 -.., 12 (2.4) 

where y- is a generalized motion variable (generalized speed or generalized coordinate) of 

St- 

Next, we imagine that the set of differential equations (4) is solved for y,‘ as 
functions of time. We note then that are the generalized motion variables for the 
variable mass system S at time t = t k . We now consider some quantity that characterizes 
the motion of ^ in some way, and is therefore obtainable from the generalized motion 
variables y*, or is simply one of them. An example of a good candidate for is the 
magnitude of the velocity of the mass center of S*. Another example is a component of the 
angular velocity of the rigid body B that is a part of S k . We will call a characteristic 
motion variable for S k . Once y * is known for all times, y k is known, and can be plotted as a 



function of time. Now, suppose we solve Eqs. (4) for various values of k, and thus 
determine ^ for several constant mass systems; and all of these ^ (* = 1, 2, .... m) are 
plotted as functions of time on the same scale, but staggered for clarity as shown in Fig. 2. 

Now, we consider the motion characteristic, v, of the variable mass system S, that 
corresponds to y_ k . If, for example, represents the speed of the mass center of S k , then v 
would represent the speed of the mass center of S, keeping in mind that different material 
particles make up S at different times. Because the dynamic behavior of S at some instant 
of time tj is, in fact, the dynamic behavior at time t } of £j, it is clear that v is given as a 
function of time by the curve labeled F in Fig. 2; or, more precisely, the projection of F 
onto the plane y.,-t (see v in Fig. 2). The task before us can be viewed as the determination 
of an efficient method for generating the differential equations whose solutions lead directly 
to the curve v of Fig. 2. The route to accomplishing this task will now be delineated. 

We consider once more Eqs. (4), which are the equations of motion for the constant 
mass system S k . First, we note that this equation has the same form for various values of 
k. The only items that change with k are system parameters such as mass and moments of 
inertia. Setting t = t k in Eqs. (4) yields a set of algebraic equations that produce y r */_* 
(r=l,2,...,12); these quantities are equal respectively to y/rj, y 2 (t,J, ...y n (t,), the time 
derivatives of the corresponding generalized motion variables for S evaluated at time t = r*. 
One of these (or some function of these) represents the slope of curve v of Fig. 2 at r = t k , 
and is plotted as point P k in Fig. 3. The above process can be repeated for The system at 
time t t , the system £ 2 at time t 2 , etc., and the relevant results are used to complete the plot 
of Fig. 3. The equations of curves such as P in Fig. 3 are the equations of motion of the 
variable mass system. Since the points of such curves are generated from the equations of 
motion of the various constant mass systems, the equations of motion of the variable mass 
system S, have exactly the same form as the equations of motion of a typical constant mass 
system. However, to apply such equations correctly to variable mass systems, care must 
be taken to interpret mass and inertia parameters correctly. At any given instant of time. 



these parameters take on their values for the corresponding constant mass system for the 
instant under consideration. 


TO)t 



Fig. 2.2 Characteristic Variables 


2.3 Dynamical Equations 


To obtain the dynamical equations of the full variable mass system, S, all we need 
to do is to derive the dynamical equations for a typical constant mass system. We will now 
do this using Kane's equations given in Eq. (1). 
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Fig. 2.3 Graph of y/t) versus time 


2.3.1 Useful Kinematical Quantities 


In order to keep the mathematical developments that are going to follow relatively 
compact, and to avoid symbol definitions that are interspersed throughout the document, 
we give and define, in Table 1 , a set of symbols that will be used repeatedly in the 
remainder of this document. 

A typical constant mass system, 2^, is shown in Fig. 4. A possible choice of 
generalized speeds, Ur, for the system is 


I co • b r (r =1,2,3) 
u r = ( , 

\\ s • b r _ 3 (r = 4, 5, 6 ) 


(2.5) 



Table 2.1 - List of Symbols 


A symbol in boldface type signifies a vector, and a boldface symbol with a tilde above it 
represents a dyadic. Several of the symbols defined here are shown in Fig. 4. 


•N -an inertial reference frame 

• B - the boundary of the variable mass system 

• R - the solid portion of the variable mass system 

• S - the variable mass system enclosed by its boundary B 

• S* - mass center of the system 

• O - an arbitrary point of R 
•P- a generic point of the system 

• v - velocity of P in an inertial reference frame 

• to - angular velocity of R in an inertial reference frame 

• a - acceleration of P in an inertial reference frame 

• v r / a, - velocity / acceleration of P relative to R 

• v° / a° - velocity / acceleration of O in an inertial frame 

• v s* / aS* _ velocity / acceleration of S* in an inertial frame 

• v?* / a?* - velocity / acceleration of S* relative to R 

• bi (i= 1,2,3) - a dextral set of mutually perpendicular unit vectors fixed in R 

• r - position vector from O to P (see Fig. 3) 

• r* - position vector from O to S* (see Fig. 3) 

• a subscript, u r , attached to any v or to indicates the corresponding partial velocity 

• p - position vector from S* to P 

• I.,- central inertia dyadic of the system 
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The velocity of a generic point, P, of the system is 


v=v° + G)xr + v r 


( 2 . 8 ) 


Because the motion of the fluid particles relative to R is assumed to be described by known 
functions of time, the partial velocities [see Kane and Levinson (1985)] of P are given by 


v — y° 
V “r V “r 


+ o»,, xr 


(2.9) 


Now, the velocity of the system mass center can be expressed as 


yS — yO + Q) x r* + 


( 2 . 10 ) 


and so, the corresponding partial velocities are 



( 2 . 11 ) 


From (11), 


XfO — yiS 

V U, - v u. 


to Ur xr’ 


( 2 . 12 ) 
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so that (9) and (12) give 

*u r mv C + t0 Ur x P (2-13) 

where p is the position vector from 5* to P. 

The acceleration of the generic point, P, of S is 

a = a° + axr + ©x(toxr) + 2o)xv r + a r (2.14) 

Similarly, 

a 5 * = a° + a x r* + co x (to x r*) + 2 <o x v£* + a?* (2. 15) 

(14) and (15) can be combined to give 

a = a s * + ax p + q)x(cdx p) + 2cox(v r - v?*) + (a r -a^*) (2.16) 

2.3.2 Generalized Inertia Forces 

The contribution of matter contained in an elementary parallelepiped located at P 
(see Fig. 4) to the generalized inertia forces of the system S* is 
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- dm a v Mr 


(2.17) 


where dm is the mass of matter inside the parallelepiped. (13), (16) and (17) yield 


(F* )p = -dm [a s * + axp + ©x(mxp) + 2cox(v r - v?*) + (a r - a?*)]- v 
- dm [a 5 * + axp + (ox(o)xp) + 2(ox(v r - v?*) + (a r - a?*)]- (o> Mr x p) 


Hence, the generalized inertia force on S* corresponding to the generalized speed, u r , is 



a s + axp + (ox((oxp) + 2 




+ axp + a>x(coxp} + 2 


(OX 


(OX 


( v r-Vr*)+( a r 
( v r- v r) + (a r 



dV 



(2.19) 


We note here that the integrals in (19) are volume integrals that are taken over the region 
enclosed by £ k . We recall that B* always contains all the material particles found inside B 
at time t k . coincides with B at the instant t k , but is different from B at any other time. 

We now return to the determination of the generalized inertia force, F* From 
Fig. 3 and the definition of mass center. 
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- r * ) dm = 0 


Similarly, 



( 2 . 20 ) 


( 2 . 21 ) 


Using (20) and (21) together with the facts that 


and 



px(ax 


p)Jd>n =1 


a 


( 2 . 22 ) 


px G>x(a>xp)p>n =J^ a> x |p x (a> x p] jdm =a»xl-a> 


)B k 


(2.23) 


(19) is reduced to 



s s 
-mv u a 



I -a + cuxlco + 2 




(2.24) 


Hence, for r = 1, 2, and 3, 
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I-a+o)xIco + 2 


L/[p x ( a,xv ^ + / B /(p xa ^] 


(2.25) 


and, for r = 4, 5, and 6, 




S 


* 


(2.26) 


2.3.3 Generalized Active Forces 

The force per unit volume, F^, acting on the generic particle P of the system can be 
written as 


F p = F£ + Ff 


(2.27) 


where F£ comes from forces external to the system, and ff is from the internal forces 
acting on P. The generalized active force, F n for the system has the form 
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Fr = 



dV 



Now, if the internal force exerted on any particle P of the system S* by another particle Q 
of Sk is assumed to act along the line connecting P and Q, then 


and 



(2.29) 


J Bt (pxFf)dV = 0 


(2.30) 


In that case, the generalized active force can be written as 


F r =v 


s* 

U r ‘ 


F + tD Ur • M 


(2.31) 


where F is the resultant external force on the system, and M is the sum of the moments of 
all the external forces on the system about the system mass center. (5) and (31) then give 
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F r = co Ur • M for r = 1, 2, 3 (2.32) 

and 

F r =v{£- F for r = 4, 5, 6 (2.33) 


2.3.4 Equations of Motion 


Recalling that Kane’s equations have the form given in (1), we now combine (25) 
and (32) to obtain the equation of attitude motion as 


I ■ a+toxIai + 2 


p |p x |cd x v r )j dV+ p(p x a r J dV = M (2.34) 


Similarly, the equation of translational motion comes from (26) and (33), and is 


ma* = F 


(2.35) 


At this point, there arc two main obstacles that prevent (34) and (35) from being 
useful as given above. First, the motion of fluid particles within was assumed known. 
This means that both the velocity and acceleration fields within jj* are known as functions 
of time. In reality, none of these functions is known, but reasonable guesses can be made 



for the velocity field within B; that is, the velocity distribution inside B* for the instant t = 
t k , when the constant mass system coincides with the variable mass system S. On the other 
hand, we have no handle over the velocity distribution of any particles outside B; and some 
such particles would normally be within B* at instants of time different from f*. The other 
problem is that we have no way of estimating the acceleration field - not even within B. 

To circumvent these problems, two important measures are taken. One is that 
attempts are made to convert all accelerations that appear explicidy in (34) and (35) to time 
derivatives of velocity. If done properly, this would take care of the second problem 
above. The other measure is that ways arc found to convert any volume integral that 
contains a velocity term, and is taken over the region Bk to a volume integral over B, where 
velocities can be estimated. This measure would resolve the first problem. We now show 
how these measures can be applied to Eq. (35). 

Substituting for a 5 * from (15), (35) becomes 


m 


a° + axr* + wx 


(ioxr*j + 2<oxy* +a* 


= F 


(2.36) 


It is convenient to express (36) as 


m 


a° + 


axr* + e>x(coxr*] +2 p (<» x v r j dV + pa r ^V = F 


(2.37) 


We then re-write the last term on the left hand side of (37) as 
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C C *dv c 

Lr'" = k <> ^ dV= «k (,v ' dV a38 

where the left superscript, R, indicates the reference frame in which the time derivative is 
taken. Next, we substitute (38) into (37) and evaluate the new expression at time t k , when 
5* is identical with S, and obtain 


m 


a° + a x r* + to x i 


R 

(" x r ')|, . , k + 2 \ B 0 (“ x v ') M + p v ' M 


=m 

t = t k (2.39) 


Because the last term on the left hand side of (39) contains a time derivative outside the 
integral sign, the limit of the integral cannot be changed to B as was done for the previous 
term. We are thus faced with the problem of taking the integral of an expression containing 
a velocity term over a region in which the velocity field is not known - the first problem 
discussed above. Fortunately, this dilemma can be resolved by means of the Reynolds 
Transport Theorem, which gives in this case. 



p \ r dV 


l/ = r* 




\t=t k 


(2.40) 


where S is the surface area of B, and n is a unit vector normal to B and pointing outwards. 
Eqs. (39) and (40) thus give 
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(2.41) 


where m is the mass of 5* st time t k , that is, the mass of S at time t k . Because (41) has the 
same form for all values of t„ the equation of translational motion for the full variable mass 

system is 


m 


a° + axr* + o>x 


a>xr +2 


J^ptcxv,)^ 


*d f 

+ dtjfl p 


v r dV + j^p v r (v r -nJdS = F 


(2.42) 


The equation of rotational motion for the variable mass system is obtained in a 
similar way from (34). First, we have that 



(2.43) 



Then we invoke Reynolds Transport Theorem once again to convert (43) to 




(2.44) 


Thus, the equation of rotational motion becomes [from (34) & (44)] 



(2.45) 


where the inertia dyadic that appears in (45) is that of matter that is within B at the instant 
under consideration. 

The vector equations of motion, (42) and (45), obtained through the application of 
Kane’s equations, are identical to those obtained by other authors [ Meirovitch (1970), 
Wang (1993)] using Newton-Euler formulation. 
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2.3.5 Another form of the Equation of Attitude Motion 

Another very useful form of the vector equation of rotational motion can be 
obtained as follows. First, we go back to the equation of attitude motion of a representative 
constant mass system as given in Eq. (34), and consider the expression 


A - K ih t [ px<<ox H = *iL t " [ p x <" x H ^ 


(2.46) 


This expression can be expanded in the following way : 

A»^J c Jpx(©xp)]dm = J^ J Jp x (< d x p)j dm 

(2.47) 

= J [v r x(mxp)]^ + J fl jpx(axp)]dw + j^Jpx(©xv r )J<im 


The last term on the right hand side of (47) can be added and subtracted from (47), then by 
using the following equality 


[v r x (co x p)] - [p x (co x v r )] = [co x (v r x p)] (2.48) 


we have 
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A = J^ x (v r x p 

On the other hand. A, as given in (47), can be evaluated at t = t k and expanded using 
Reynolds Transport Theorem to yield 


dm+ I a + 2 px(a>xv r 


(2.49) 



(2.51) 
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Eqs. (34) and (51) lead to the following alternate form of the vector equation of attitude 
motion: 


a+toxi(o+ — a> + J^p|p x(<Dxp)](v r -n)dlS 



xv r)^ + J 5 P( p x v r ) (v r n) <iS = M 


(2.52) 
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CHAPTER 3 

EFFECTS OF MASS VARIATION 


The equations of motion derived in chapter 2 [see (2.42), (2.45) and (2.52)] are 
quite complex; and very little progress can be made with these equations without some 
simplifying assumptions to reduce their complexity. In this chapter, these equations are re- 
examined. The significance of each of the terms is evaluated, and , the circumstances 
under which some of these terms can be dropped are explored. The goal is to arrive at a set 
of equations that capture the salient features of the system under study, and that are 
nevertheless simple enough to render further analysis tractable. 


3.1 Translational Motion 

The vector equation of translational motion [see (2.42)] is reproduced here as 


m 


a° + axr* + tDx 


(a> x r*) + 2 J^p(® xv r) 


dV 


+ SL fir '‘ iv+ Is fi, '( r ' n ) ds=F 


(3.1) 



This equation can be re-written as 


m 


a° + a x r* + a> 


x [a> x r*) 


— F + F r +Fi + Ft 


(3.2) 


where 



(3.3) 


Fl 



py r dV 


(3.4) 


and 



(3.5) 


If v r is identically zero everywhere, then F c = F L = F r = 0, and we recover the equations 
of translational motion for a rigid body. Thus, mass variation appears to augment the 
“external forces” on an “equivalent rigid body” by three terms. Fc is often referred to in 
the literature as the Coriolis force, since it derives from the Coriolis component of the 
acceleration. Fl is the rate at which the system’s linear momentum relative to B decreases 
with time because of particle motion inside B. Ft represents the rate at which relative linear 
momentum is lost across the boundary B, and is often referred to as the thrust vector in 
rocket applications. 

If those particles of the system that can move relative to B are allowed to move 
within B but do not cross the boundary B, then F T becomes zero. F x would be non-zero 
but negligible if either a very small percentage of the system’s particles is allowed to cross 



B, or those particles that cross B do so at very slow rate. In other words, the thrust vector 
can be neglected whenever the amount of matter that is lost or gained per unit time is small. 
What matters is the rate, not the total amount of matter lost or gained. 

Note that we can use Reynold’s transport theorem to expand the Coriolis force into 


F c = - 2 to x J py r dV=-2tox-^J prdV 


=-2 <D>d 


"4 /, » r + 1 P r < V ' "> ^ = F ‘l + F «2 


(3.6) 


where 


^. = - 2 “ x s// rdV 

and 

F C2 =-2 to x £ p r (v r n) dS 


(3.7) 


(3.8) 


There are situations where v r can be considered negligible inside the boundary B 
but not at an exit from B. An example is an inflated balloon with a hole. Gas motion 
inside the balloon is hardly noticeable while gas exit velocity at the hole is substantial. In 
cases like this, F L and the first term of F c (F c/ )are negligible, but ¥ C2 as well as the thrust 
vector survive. As a matter of fact, this is not an unreasonable assumption for practical 
systems such as rockets, and can provide a way of rendering the equations of such systems 
tractable, since the details of internal gas flow can be neglected. Even if v r is not negligible 
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within B, but the particles within B can attain some type of steady state in their motion 
relative to B, then F L and F CJ are once more negligible. After ignition, rocket systems 
quickly attain an approximate steady state, and so, F t and F c , can be ignored for such 
systems. 

In summary, the vector equation of translational motion for variable mass systems 
can, in many cases be reduced to 


m 


a° + axr* + o>x 


(a> x r*) 


— F + Ft +F ( 


Ci 


(3.9) 


where F is the external force, F r is the thrust vector, and F a is part of the Coriolis force 
and is given by (8). F c and F r can be dropped completely only when the rate of mass loss 
is negligible. 

3.2 Rotational Motion 

Two forms of the vector equation of attitude motion were developed in chapter 2 
[(2.45) and (2.52)], and are reproduced below. 


I ■ a+ ©x I •© + 2 




<JV- 


+ *' Vr ) dV + l P ( P XV r)( V r")^ = M 


(3.10) 
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I - a+coxI-co + ^-^J © + J^p|px(cDXp)J(v r -n)dS + 

R 

J fl p|©x(px v r )] dV+ -£f s p (p x v r ) dV + p (p x v r ) (v r n) dS = M 


(3.11) 


Both of these equations reduce to rigid body rotational equations if v r = 0 within B as well 
as on B. Thus (10) can be written as 


I • a + © x I co = M + M c + M h + M x 


( 3 . 12 ) 


where 


M c = -2 



(3.13) 


is the Coriolis moment. 


M w = - 



dV 


(3.14) 


represents the rate of decrease of the system’s relative angular momentum inside B, and 
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M r = - j s P[v x v r J(v r -nJdS (3.15) 

is the rate of loss of relative angular momentum across the boundary, and is also equal to 
the moment of the thrust vector about the system mass center. 

In those situations where it is reasonable to assume that the motion of the fluid 
phase relative to the solid phase has axial symmetry, and when there is no whirling motion, 
we have that 


p(px\ r )dV = p(px\ r )dV=Q 
jBlr JB 

k (3.16) 

This is so because, for every particle P of the system with position vector p and relative 
velocity v n there exists another particle P’ of position vector p’ and relative velocity v’ r 
such that the vectors pxv r and p’xv’ r have the same magnitude but opposite directions. 
And, since axisymmetry also implies that the immediate neighborhoods of P and P’ have 
the same mass density, Eq. (3.16) follows. This equation immediately leads to 

M h = M t = 0 (3.17) 

We arrive at the same conclusion for systems where v r can be considered negligible inside 
the boundary B but not at an exit from B, such as rockets. M H is also negligible for steady 
state relative motion of the fluid particles. 


From (10) and (1 1) the Coriolis moment can be expanded to 



(3.18) 


M c = - 



>)](v r -n)ttf + J^pJ«>x(p 



dV\ 


The last term on the right hand side of (18) vanishes for axisymmetric motion as well as for 
negligible internal flow. The second term on the right hand side of (18) is often referred to 
as the jet damping moment, because it has been shown [see Dryer (1963)] to have 
attenuating effects on the angular rates in some types of rocket systems. To evaluate this 
surface integral, it is necessary to know the geometric shape of those parts of the boundary 
where particles are allowed to exit or enter the system, as well as the velocity profile at 
these locations. This moment is thus very much dependent on the system’s geometry. In 
the case of rockets, for example, the longitudinal dimension of the combustion chamber 
(for solid rockets) as well as the exit nozzle radius have much to do with the impact of this 
term on the system’s attitude motion. 

The first term on the right hand side of (18) captures the contribution of inertia 
variation. We thus conclude that jet damping moment and the moment due to the changes 
in inertia properties have dominant effects on the attitude dynamics of the system. The 
relative importance and the interaction between these two moments essentially determine the 
character of the attitude dynamics of variable mass systems. There is advantage in using 
the reduced form of (11) for the study of rotational motion. If the assumptions of 
symmetric and/or negligible internal motion are made, the most “troublesome” terms will 
have dropped out. The only term in this equation that would contain v r is a surface integral 
over B. v r -n is zero everywhere on the surface of B except at those places where fluid 
particles can enter or leave the region delimited by B (the nozzle exit plane in the case of 



rockets). In general, v r n can be approximated relatively well at these locations; hence the 
surface integral can in fact be evaluated in closed form. 

In the study of the dynamic behavior of variable mass mechanical systems, the most 
important forces appear to be the thrust vector and the Coriolis force. The jet damping 
moment and the moment due to inertia variation are the dominant moments. These 
quantities should be included in any dynamic studies of variable mass systems to ensure 
that meaningful predictions can be made from the analysis. 



CHAPTER 4 


ATTITUDE MOTIONS 
OF A VARIABLE MASS CYLINDER 


4.1 Introduction 

In this chapter, we begin to narrow the scope of our study to an important and 
useful class of variable mass systems - space rockets. We start this process with an in- 
depth study of the attitude motions of a variable mass system that is initially a solid right 
circular cylinder, which loses its mass continuously through combustion as it moves 
around in space. The cylinder problem is a very important and useful problem for several 
reasons. First, it models rocket-type systems in a simple enough way that the equations of 
motion become relatively tractable. Generally, space vehicles are designed to be more or 
less axisymmetric is shape, and, a cylinder is actually a good - albeit rough - 
approximation for such systems. Hence, a thorough study of the variable mass cylinder 
problem is a useful exercise in that it provides a means of performing tractable analytical 
studies of rocket-type systems, and can lead to great insight into the dynamic behavior of 
this class of variable mass systems. 



The idea of using a cylinder to study the attitude behavior of variable mass systems 
was originated by Eke and Wang (1995). In their work, they solved completely in closed 
form, the equations of motion of a variable mass cylinder for several bum scenarios. They 
then extracted useful qualitative as well as quantitative information about the attitude 
behavior of the cylinder from the analytical expressions of the solutions of the equations of 
motion. The approach here is different The strategy here is to develop an analytical 
method that yields a wealth of fundamental information about the behavior of the system 
without actually solving the equations of motion. The main motivation for this approach is 
that such a strategy can be applied to more complicated models of variable mass systems 
for which closed form solutions cannot be found. 


4.2 Attitude Equations for the Cylinder 

Consider a variable mass system of the rocket type that is initially a right circular 
cylinder as shown in dotted lines in Fig. 1. This system is given an initial angular velocity 
G\» and is then allowed to move freely in a torque-free environment as it loses its mass 

continuously through combustion or similar processes. The intent is io perform qualitative 
and quantitative studies of the attitude motions of such a system. To begin, we define as 
control region, the space enclosed in dotted lines in Fig. 1. Matter within this region at any 
given instant is considered part of the system at that instant 

The attitude motions of such a system are governed by Eq. (2.52), which can be 
simplified, as discussed in chapter 3, by making the following assumptions: (a) the 
burning process proceeds in such a way that the solid unbumed part of the cylinder remains 
symmetrical about the original cylinder axis at all times; (b) the fluid products of 
combustion move in an axisymmetric manner relative to the solid unbumed part of the 



cylinder, (c) whirling motion of the products of combustion relative to the solid portion of 
the system is negligible. With these assumptions, (2.52) reduces to 
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Fig. 4.1 Variable Mass Cylinder 


The axisymmetry assumption allows us to express the instantaneous central inertia dyadic 
of the system as 


I-/(b l l>i + b2b2j +/ z b3b3 


(4.2) 




where b/, b 2 , bj are unit vectors fixed in B and directed as shown in Fig. 1; and 7 and 7, 
are respectively the transverse and axial central inertia scalars of the system. If the inertial 
angular velocity of B at some general instant of time is taken to be 

o> = n^b ! + dfyb 2 + (43 


then, the angular acceleration is 


a= a^bj + a^b 2 + t^b 3 


(4.4) 


and the various terms of Eq. (1) become 

I • a = 7 + 0^)2) + / z ^3 


(4.5) 


cdxI • a>= (7 z -/)<q (fi^bi-fi^b^ 


(4.6) 



a>= 7 (tqjb ! + ojjb^) + 7 Z (ap 3 


(4.7) 


In order to evaluate the last term of (1), an assumption must be made concerning the 
velocity field of the fluid particles of the system at the surface of the control region. Here, 
we assume that the vector v r is zero everywhere on B except at one of the ends of the 
cylinder — the right end, say. And there, 

V" = “('■> (4.8) 


because of our axisymmetry assumption. If it is further assumed that 



u(r) = u = constant 


(4.9) 


then 


f. 


p(px(<Bxp)(v r -n)] dS = -ni (( 2 +^R 2 )(Ciip l + ^ 2 ) + ^ 2 ^ b 3 


(4.10) 


where m is the instantaneous mass of the system, and /(f)is the distance of the mass center 
from the right face of the shell, B. By substituting (5), (6), (7) and (10) into (1) and 
setting M to zero, we obtain the scalar equations that describe the torque-free attitude 
motion of the variable mass cylinder, and they are 


/4 + (/ 2 -/)^^ + |/-m|/ + /? 2 /4j|ai = 0 
/ 2 -»- { / z - /? 2 /2 


(4.11) 

(4.12) 

(4.13) 


4.3 Attitude Stability 

Multiplying (1 1) by co*, (12) by d) y and adding the two resulting equations leads to 


1 -d- (a)* + cd£) + [ i - m (t 2 + R 2 / 4 ) ] (cal + = 0 

4W UJL 


(4.14) 
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We now let 


0 = i r - m R 2 / 2 (4.15) 

<J> = I - rh (/ 2 + R 2 / 4 ) (4.16) 


and 


(Dt = G& + (£$ 


(13) and (14) can then be written more compactly as 


j d(0 2 

iz ~dT 


+ 0 co z = 0 


(4.17) 


(4.18) 


and 


I ^ + 2 <]> to, = 0 
at 


(4.19) 


If the functions 0{t) and <|)(t) are known, (18) and (19) can serve as basis for the 
determination of the essential features of the rotational motion of the variable mass cylinder. 
We note that these equations arc uncoupled, so that the spin rate has no effect on the 
transverse angular rates. We also note that if 0(t) = 0, then the spin rate remains constant at 
its initial value throughout Similarly, if <|)(t) = 0, the magnitude of the transverse angular 
velocity does not change during the bum. 

Equations (18) and (19) lead to 

to* = CDz 0 exp 



(4.20) 
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or 


0), = <Df o exp 




CDjy — CDjy o exp 




(4.21) 


(4.22) 


where (Oxy = Va^ is the magnitude of the transverse angular velocity, and the subscript 0 
indicates initial value; that is, value at ignition. We shall use the subscript <» to indicate 
values at burnout. 

Looking back at (15) and (16), we observe that each of the functions 0 and <|> is 
made up of two parts. For example. 


6 = 01 + $2 


(4.23) 


where 


= 1 


(4.24) 


and 


te = -m(f 2 + R 2 /4) (4.25) 

Because we are concerned here with a situation where mass is being lost, both i and rii 
are negative. Hence, <f>i remains negative throughout the bum while 62 is positive 
throughout. Thus, <j> can be positive, negative or zero. <|>i is contributed by the change in 



system transverse inertia, and because it is always negative, it will tend to cause a ) xy to 
diverge [see (22)]. On the other hand, <(>2 is, in this case, the so-called jet damping term, 
and it does in fact attenuate co x> by virtue of the fact that it is always positive. Although <]>i 
is determined by the bum geometry of the cylinder, <J >2 depends on the size of the cylinder 
as well as the velocity distribution across its exit plane. All the statements made above 
about <Kt) apply equally to 0(t) since (15) and (16) have the same form, and (20) is similar 
to (22). 

It is possible to extract important qualitative information about the attitude motions 
of the system from a detailed examination of the functions <ji(t) and 0(t). From (20) and 
(22), it is clear that for 0(t) and <f)(t) > 0, both co* and (O xy approach zero from any initial 
conditions. On the other hand, for 0(t) and <^t) < 0, both to* and (0 xy diverge. Profound 
changes in angular velocity can occur when variations are made in 0(t) and <b(t). 

To make further progress with our study of the attitude behavior of the variable 
mass cylinder, we will now proceed with a close study of the functions 0(t) and <|>(t). 

These functions require that the system mass, inertia scalars and some geometric properties 
be known as functions of time; and these parameters can be determined if the bum scenario 
is known. We will therefore consider, as did Eke and Wang(1995), four bum scenarios : 
the uniform bum, the end bum, the radial or centrifugal bum, and the centripetal or anti- 
radial bum. In uniform bum, the assumption is that the cylinder is ignited simultaneously 
everywhere inside it at time to- It then bums in the same uniform manner everywhere, so 
that the external dimensions remain the same throughout, but the density of matter within 
the cylinder decreases uniformly in the same manner everywhere within the cylinder. The 
density is thus the same function of time everywhere inside the cylinder, and the products 
of combustion are expelled at one end of the cylinder. An end burning cylinder bums from 
one of its ends to the other, in such a way that the intermediate shape of the system is 
always a cylinder of the same radius but decreasing length. Once more, the products of 
combustion are expelled from one of the ends - the burning end. In radial bum, the axis of 
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the cylinder is ignited at time to, and the system then bums radially outwards in such a way 
that the intermediate shape is a cylindrical pipe. Counter-Radial or centripetal bum is the 
reverse of the radial bum. The cylinder is ignited at its periphery (but not the ends) and 
bums radially inwards, with the intermediate shape being a cylinder of constant length but 
of decreasing radius. 

4.3.1 The Uniform Burning Cylinder 



We now consider the case of a cylinder in uniform bum with uniform velocity 
profile at the exit plane, as shown in Fig. 2 above. With the dimensions shown in the 
figure, the inertia scalars are, in this case, 


I = m ( R 2 /4 + h 2 /3 ) 


(4.26) 


and 
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Hence 


I z = m R 2 /2 


(4.27) 


I = m ( R 2 /4 + h 2 /3 ) 


(4.28) 


and 

iz = m R 2 /2 (4.29) 

From (15) and (29), we have that 

8(0 = 0 (4.30) 

at all times, and so, the spin rate, Gfe, will remain constant at its initial value. 

From (16) and (28), and recalling that l in (16) is the same as h in Fig. 2 for this case, 

<|>(0 = -|mh 2 >0 (4.31) 

Because <t>(t) > 0 at all times during the bum, the magnitude of the transverse angular 
velocity, C 0 x y , decreases exponentially to zero. 

4.3.2 The End Burning Cylinder 


Fig. 3 shows the dimensions and intermediate shape of a cylinder in end bum. The 
auxiliary spatial coordinate z is employed to characterize unbumed propellant inertia 
properties, including mass, moments of inertia and the location of center of mass. The 
mass and inertia properties and their time derivatives are 
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m = 2 p re R 2 z 

(4.32) 

m = 2 p 7t R 2 z 

(4.33) 

I = m ( R 2 /4 + ztf'i ) 

(4.34) 

i = m ( R 2 /4 + z 2 ) 

(4.35) 

I z = m R 2 /2 

(4.36) 

i z = m R 2 /2 

(4.37) 


h 



2h-z 


Fig. 4.3 Cylinder in End Bum 


Furthermore, 



50 


f = ( 2h — 2z) + z = 2h — z 


(4.38) 


so, from (16) and (38), 


<]>(t) = m ( R 2 /4 + z 2 )-m[(2h-zj 2 + R 2 A ] 
= mh 2 [4(z/h)-4] 


(4.39) 


At ignition t = to and z = h, so that 


<t>0 = 0 


(4.40) 


At burnout t = too and z = 0; hence 


4>oo = - 4 m h 2 > 0 


(4.41) 


From the mass continuity equation, we also have 


m = -pnR 2 u<0 


(4.42) 


Based on the assumption of uniform velocity profile at the exit plane, we conclude that m is 
constant during the time interval [to, to©]. We also know that z is negative in this interval. 
From (39), we have that 


d^/dt = 4 h m z > 0 


(4.43) 
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between ignition and burnout. Thus, <)>(t) increases monotonically from 0 to the positive 
value given in (41). We have a situation where <J>(t) has the form shown in Fig. 4. 



Fig. 4.4 The shape of <J)(t) for End Bum 


We conclude from (40), (41) and (43) that the transverse angular rate, G) X y, stays 
approximately constant during the initial phase of the bum, then decreases exponentially to 
zero. 

As for the spin rate, (15) and (37) indicate that 0(t) = 0 for all t, so that the spin rate 
is again constant 

Having established that the transverse angular velocity of an end burning cylinder is 
always damped, our next task is to examine the effect of the cylinder’s initial geometric 
configuration on the lateral response. To accomplish this task, let us define 

'F = / I(t) (4.44) 

It is clear from (22) that this factor actually governs the speed with which the transverse 
angular speed approaches zero from any initial value. The larger the magnitude of V F, the 
higher the rate of decay of lateral motion. 

For end bum. 
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y_ < t>_rii (4hz-4h 2 ) 
1 m (r 2 /4 + Z 2/3) 


(4.45) 


By virtue of (32) and (42), (45) can be rewritten as 


2u(l-z/ h ) 

z[l/4(R/hf+l/3(z/h) 2 ] 


(4.46) 


To compare the behavior of cylinders of the same length but different initial external radii, 
we allow R to vary while keeping z, h and u invariant. A close inspection of (46) indicates 
that increasing R / h results in smaller values for V F. Hence, the transverse angular velocity 
magnitude, G) xy , converges more slowly as the ratio R / h is increased. In other words, the 
transverse motion of prolate cylinders in end bum is highly damped while that of oblate 
cylinder is only very slightly damped. In the limiting case, as R/ h — » «>, = 0, the 

motion is totally undamped as in the torque-free motion of a rigid cylinder. We conclude 
then that the amplitude of lateral oscillation of a variable mass cylinder in end bum can be 
influenced by the choice of its initial shape. Fig. 5 is obtained by numerical integration of 
(22) and demonstrates all the above properties. 


4.3.3 Radial Burn 

We recall that in radial or centrifugal bum, the cylinder starts burning on its axis, 
and the bum propagates radially outwards with time. The intermediate configuration of 
unbumed propellant is a hollow cylinder, as shown in Fig. 6, and the mass and inertia 
properties in terms of auxiliary spatial coordinate r have the form 




54 


m = 2phn(R 2 -r 2 ) 

(4.47) 

I = m [l/4 (r 2 + r 2 ) + h 2 /3] 

(4.48) 

+ 

'Si 

II 

(4.49) 

Hence 


m =- 4 p h Jt rr 

(4.50) 

It is relatively easy to show that 


i = m ( x^2 + h 2 /3 ) 

(4.51) 

and 


i 7 = m r 2 

(4.52) 

To look at the behavior of the spin rate, we examine 


e(t) = i z -(m/2)R 2 = mR 2 [(r/Rj2_ 1 12 ] 

(4.53) 


Assuming that at time t = 0, r = 0, we have 


6(t)lt = 0 = -mR 2 /2>0 


(4.54) 



and at the end of the bum ( t — » too, r = R), 


e(t)|t-Moo = mR 2/ 2<0 


(4.55) 


From (53), we have that 


d0/dt = 2mrr<O (4.56) 

between ignition and burnout. Because r varies from 0 to R in the time interval [ 0, too] and 
from the physics of the problem, it is clear that r > 0 in the given interval. 

We thus have a situation where 9(t) has the form shown in Fig. 7 

act) 



Fig. 4.7 Shape of 0(t) for Radial Bum 


The switchover point, t*, at which the sign of 0(t) changes from positive to negative is 
obtained by setting 0(t) = 0 in (53) : 


r /R = yT7T = 0.707 


(4.57) 
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This implies that the spin rate decreases in the first stage of the bum, when r / r <; 0.707, 
but increases exponentially in the later phase of the bum when r/R > 0.707; a fact that is 
clearly supported by the plot shown in Fig. 8. This plot is obtained from numerical 
integration of (20) for centrifugal bum. 



Fig. 4.8 Spin Rate Curve for Centrifugal Bum 
We now discuss the evolution of the magnitude of the transverse angular velocity, 

(O xy . 


<j)(t) = i- m(f 2 + R 2 /4) = mR 2 [(1/ 2) (r / r)?- ( 2/ 3)(h / RJ 2 -1/4] (4.58) 


Hence 


<Kt)|t=o = 4>o> 0 


(4.59) 



and 


^ = ihR 2 [l/4-(2/3)(h/Rf] (4.60) 


From (58), 


d<j>/dt = mrr<0 (4.61) 

The sign of <t>« depends on what can be described as the cylinder’s “shape factor”- the ratio 
of die initial value of its diameter to its length. Hence, whether is bounded or 
unbounded depends on this shape factor. To determine the value of R / h that separates die 
“stable” region of ct> xy from the “unstable" region, we set <K0 to zero in (60), and obtain 

R / h = V8TT = 1.63 (4.62) 

Relevant “pictures” for <j>(t) are shown in Fig. 9. 



Fig. 4.9 Shape of <(>(t) for Radial Bum 
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Hence, the magnitude of the transverse angular velocity of a cylinder in radial bum always 
decreases with time in the early phases of the bum. If the ratio of the unbumed cylinder’s 
diameter to its length (shape factor) is less than that of 1.63, the decrease of o) Z)t continues 
all the way to burnout. On the other hand, if the cylinder’s shape factor is greater than 
1.63, then there comes a time during the bum when © xy levels off, and then starts 
increasing exponentially through burnout. The initial shape of die cylinder thus has a 
critical effect on the cylinder’s lateral attitude motion for radial bum. Pencil-shaped 
cylinders will tend to be stable in radial bum while hamburger-shaped cylinders will tend to 
be unstable. We note here that the same shape factor that determines stability or instability 
in radial bum was shown to influence the rate of convergence of to X y for end-burning 

cylinders. 

It would be useful to determine the onset of instability of a> xy for cases where R / h 
> 1.63. This can be done by determining the value of r/R at which the sign of <t>(t) 
changes from positive to negative in (58). We thus have, from (58), that 

<4 - 63) 


So, the onset of instability for a) Z) , depends again on the shape factor. Once a value is 
specified for the shape factor (R / h), the value of r / R at the beginning of instability can be 
determined. For example, for 



0.913 

0.744 

0.707 


(4.64) 



We conclude, from (64), that a> xy begins to diverge much earlier for cylinders with higher 
values of the shape factor (fat and short cylinders) than for cylinders with lower values of 
R / h. Numerically obtained plots of the transverse angular velocity are shown in Fig. 10, 
and are totally consistent with the above inferences. Furthermore, the above analyses 
capture a great deal more of the vital characteristics of radial bum than do these plots and all 
previous work on the cylinder. 



Fig. 4. 10 Transverse Angular Rate for Radial Bum 




4.3.4 


Centripetal Burn 


In centripetal bum, the cylinder bums radially inwards without changing its length. 
The intermediate shape is a solid cylinder of diminishing radius, as shown in Fig. 1 1. In 
this case, we have that 


m = 2 p h 7t r 2 

(4.65) 

I 7 = m r 2 / 2 

(4.66) 


and 



(4.67) 


Fig. 4. 1 1 Cylinder in Centripetal Bum 
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The time derivatives are 


m = 4phrcrr 


U 


= m r 2 


and 


i = m ( r 2/2 + h 2 /3 ) 


For a study of the spin rate, we have 


0(t) = i z -(l / 2) (m R 2 ) = m (r 2 - R 2 / 2 ) 


00 = m R 2 / 2 < 0 


0oo = - m ( R 2 / 2 ) > 0 


(4.68) 

(4.69) 

(4.70) 

(4.71) 

(4.72) 

(4.73) 


and from (71), 


d©/dt = 2 m r r > 0 (4.74) 

since r < 0 in this case. We thus have a situation where 0(t)has the form shown in Fig. 
12 
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Fig. 4.12 Shape of 0(t) for Centripetal Bum 


Hence the spin rate increases initially, attains a peak value then decreases rapidly with the 
bum. The time of attaining the peak value of the spin rate corresponds to 0(t) = 0. That is 


r/R = 1 HI = 0.707 (4.75) 


This trend is validated by the plot of Fig. 13 below. 



Fig. 4.13 Spin Rate of Centripetal Bum [Numerical Solution] 



To assess the evolution of the transverse angular velocity magnitude, we find 


<Kt) = 1 - m ( l 1 + R 2 / 4 ) = m R 2 [( 1 / 2) (r / R j*- (2/ 3) (h / R j 2 - 1 / 4] (4.76) 


which is exactly the same as (58) for radial bum. 


<j >0 = m R 2 [l /4-(2/3)(h/Rp] 

(4.77) 

<|)« = -m R 2 [l / 4 + (2 / 3) (h/R) 2 ] > 0 

(4.78) 


and from (76) 



(4.79) 


From (77), we determine that 4>o = 0 when R / h = 1.63, <t>o > 0 for R / h < 1 .63, and 
<|>o < 0 for R / h > 1.63. We conclude from all of this that (a) as long as the ratio of the 
cylinder’s initial diameter to its length is less than 1.63, the transverse angular velocity 
magnitude decreases as the bum progress; (b) if the ratio R / h is greater than 1.63, the 
magnitude of the transverse angular velocity increases initially until it reaches a peak value, 
then it decreases with time for the reminder of the bum. Thus <o xjr never really becomes 
unbounded for centripetal bum. These facts are demonstrated below in Fig. 14. 

Our findings so far indicate that the stability of rotational motion of a variable mass 
cylinder depends on two factors : the bum pattern and the cylinder’s shape factor. Both the 
spin rate and transverse angular velocity are bounded for all shape factors when the 
cylinder is subjected to uniform, end, or centripetal bum. For radial or centrifugal bum, the 
spin rate is always unbounded, while the transverse angular speed is bounded for the shape 



factor (R / h) that arc less then 1.63, and unbounded for shape factors that are greater than 
1.63. 

All the results obtained in this section are totally consistent with, and augment 
considerably, earlier results by Eke and Wang(1995). 



Fig. 4.14 Transverse Angular Rate of Centripetal Bum [Numerical Results] 

4.4 Time Dependence of Motion and Inertia Properties 

Up to this point, all variables and parameters of interest in this study have been 
defined or determined as functions of the auxiliary variable z or r. This approach turned 
out to be a very astute and efficient analytical strategy. Through this choice, we have been 
able to avoid some of the myriad problems encountered by previous investigators of the 
dynamics of variable mass systems. For example, practically all previous investigators 
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worked directly with time as independent variable. Hence, to make progress with the 
equations of motion [ for example (1 1)-(13)], it was necessary to specify the time 
dependence of the mass (m), the location of center of mass (£) and moments of inertia ( 

I, I z ). These investigators simply made assumptions such as linear mass variation - 
m = mo - a t, and made “reasonable” guesses for the coefficients. They made similar 
guesses for l, I and I z . The problem is that these inertia quantities are interrelated; a fact 
that is generally not correctly reflected in the guesses for the independent time functions 
assigned to these quantities. All such difficulty is avoided by the clean analytical scheme of 
judiciously introducing the auxiliary variables z or r, as done in this study. 

Although avoiding the explicit use of time in our analysis so far has been quite 
beneficial from analytical point of view, it is in fact still desirable to know how the angular 
rates and even the inertia properties vary with time. For example, this enables us to 
compare analytically predicted trends with experimental results, obtained through on-board 
instrumentation. In this section, we show how some of the results obtained so far can be 
converted to time functions. 

4.4.1 End Burning Case 

We now reexamine some of the results obtained for the end burning cylinder. 

From (33) and (42), we have that 


z =-u/2 


(4.80) 


or 
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w* 

(4.81) 

which leads to 


z = - (u / 2) t + h 

(4.82) 

The inertia properties can now be converted to time functions simply by substituting (82) 
into their expressions. Thus, from (32) and (82) 

m = 2p7iR 2 z = 2p7iR 2 (h-^utj = mo + mt 

(4.83) 

where 


m 0 = 2 p ji R 2 h 

(4.84) 

and m is as given in (42). Similarly, 


I = m (r 2 / 4 + z 2 / 3) = (m 0 + m t)[R 2 / 4 + (h-ut/2)l 2 /3] 
= Io + Pt-Qt 2 + St 3 

(4.85) 

where 


£ 

ii 

(4.86) 

•6 

II 

(4.87) 


Q = i- mhu 


(4.88) 



The value of i. e. the burnout time, is obtained by setting z = 0 in (82) : 


t 00 = 2h/u 


(4.90) 


which is independent of the radius of the cylinder. The instantaneous transverse central 
radius of gyration, k, is given by 


So, from (85) 


I = m k 2 


(4.91) 


l. 2 _ R 2 . h 2 h u t + u 2 t 2 
4 1 1 12 


(4.92) 


Since 




(4.93) 


k2 = kJ-lUl-< + ^< 2 


(4.94) 


For the time interval 0 <. t <, 2 h / u, k 2 is a decreasing quadratic function of time. 


At burnout 
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k* = R 2 /4 


(4.95) 


The axial moment of inertia is 


I z = |raR 2 = ^R J (m c + m t) (4.96) 

and the axial radius of gyration is a constant. Since m and u are constants, we see that the 
transverse moment of inertia is a cubic function of time while both the axial inertia and the 
mass are linear functions of time for the end bum assumption. 


4.4.2 Centrifugal Burn 


For centrifugal bum, we have, from (47), that 


m = - 2 p h 7t = - p n R 2 u 


(4.97) 


Hence 


d (r 2 ) _ r 2 , 
dt 2 h 


(4.98) 


Integrating both sides over appropriate limits, we have 


r 2 - R 2 u 
7 h 


t 


(4.99) 
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Hence, 


where 


m = 2phjc(R 2 -r 2 ) = 2ph7t [r 2 - (r 2 u t)^(2 h)] = m G + m t (4. 100) 


= Io + B t + A t 2 



(4.101) 


Io = m 0 ( R 2 /4 + h 2 /3 ) 


(4.102) 


A = (m R 2 u/(8 h) 


(4.103) 


B= m h 2 /3 


(4.104) 


Since both A and B are negative, it is clear from (101) that I is a decreasing quadratic 
function of time. Note also that 


i = B + 2At<0 


(4.105) 


From (49) 
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h = i- m (R 2 + r 2 ) = i-(m 0 + m t)(R 2 + BluJ.) 
= Izo + 2 A t 2 


and 


1 7 = 4 A t < 0 


where 


Izo - 2 m o ^ 

Thus, I z is also a decreasing quadratic function of time. 

It is straightforward to show that the radii of gyration arc given by 


k 2 = k§ + 1 


and 


k?=k “ + ^h M 

where 




(4.106) 


(4.107) 


(4.108) 


(4.109) 


(4.110) 


(4.111) 


(4.112) 


So, k 2 and k 2 .are both linearly increasing functions of time. 
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Now, the value of burnout time U is obtained by setting r = 

R in (99). Hence 

2h 


t °°“ u 

(4.113) 

and 



(4.114) 

hi 

H 

(4.115) 

4.4.3 Centripetal Burn 


In the case of centripetal bum, similar analyses lead to 


r 2 = R 2_R2jL t 

(4.116) 

7. h 


2h 


u 

(4.117) 


m =mo + rht 


(4.118) 


I^Io + Ct-At 2 (4.119) 

with 

C = rh(^ + 1^) (4.120) 
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i=C-2 At 

(4.121) 


Iz=Izo + Dt-2 At 2 

(4.122) 

where 




D=m R 2 

(4.123) 

Hence, 

, both I and I z are decreasing quadratic functions of time. 



The radii of gyration can also be shown to be given by 



A 

i 

<3? 

n 

(4.124) 

and 




k? = k£,-^-t 

(4.125) 

Both k 2 and k 2 are decreasing linear functions of time. We also have that 



k 2 =h 2 /3 

(4.126) 

and 




O 

It 

(4.127) 


In conclusion, we observe that for the three cases examined, the cylinder mass is a 
linear function of time and the bum duration is the same for all cases. For the end burning 
cylinder, I is a cubic function of time, while k 2 is a quadratic function of time; for 
centrifugal and centripetal bums, I is a quadratic while k 2 is linear in time. All of the inertia 



parameters are decreasing functions of time except for centrifugal bum where the radius of 
gyration increases with time. 



CHAPTER 5 


MOTION OF AXISYMMETRIC BODIES WITH 

MASS LOSS 


In the previous chapter, the important problem of the attitude stability of rocket 
systems is initiated, with such systems idealized by a right circular cylinder. In particular, 
the effects of mass loss for various burning patterns are emphasized. 

In this chapter, another step in the study of the attitude motions of rockets is taken 
by considering a more sophisticated model as shown in Fig. 1 below. The system 
investigated here differs from the cylinder model in two basic ways : (i) the system is 
modeled as a general axisymmetric system ( not just a cylinder), and includes a cylindrical 
solid propellant grain as well as a constant mass, non-propellant portion; (ii) a nozzle 
arrangement is also included. The new system can be reduced to the cylinder model by 
setting the mass of the invariant part to zero. 

The objectives in this chapter are similar to those of the previous chapter, we are 
interested in the angular rates of the system, and how they are influenced by mass loss. 
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U 




Fig. 5.1 Axisymmetric Solid Rocket System 

5.1 Governing Equations 

The system under consideration in this chapter is shown in Fig. 1, and consists of 
two main parts, S and F. F represents the solid propellant of a typical rocket system. It is 
assumed to be rigid and cylindrical. S is the main part of the system excluding fuel and 
includes the payload. S is assumed to be rigid and axisymmetric, and maintains the same 
mass throughout Both S and F have the same axis of symmetry, so that the mass canters 
S* of S, F* of F and B* of the combined system, all lie on the common axis of symmetry 
z. The system loses mass continuously through a planar portion of its surface that we shall 
refer to as the nozzle exit plane, and that is located at the right end of the system in Fig. 1. 
At any given instant, the overall system will consist of everything that is bounded by the 
external boundary of S and the nozzle exit plane. 

System dimensions are characterized by the symbols given in Fig. 1. Ms 
represents the mass of the invariant part, S, and mf is the instantaneous mass of the fuel, F, 
The origin of the coordinate axes x, y, z is located at, and moves with, the composite mass 
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center B*, such that the z-axis is always aligned with the symmetry axis, and the x and y 
axes have no rotation relative to S. The solid propellant, F, is assumed to maintain an 
axisymmetric configuration throughout its bum, so that the central inertia matrices for S 
and F are diagonal. Externally applied moments, thrust misalignment effects, aerodynamic 
and gravitational moments arc assumed to be of relatively minor importance during 
burning, and are thus neglected. The system is thus assumed to move in a torque-free 
environment. The motion of the particles of the products of combustion as they cross the 
exit plane is taken to be such that the velocity distribution perpendicular to the plane is 
uniform. If it is further assumed that the motion of gas particles inside the combustion 
chamber is symmetric relative to the z-axis, and that whirling motion of the gas particles 
wi thin the system is negligible, the vector equation of attitude motion for the system is 
reduced from (2.52) to 

p[p x (<b x p)] (v r n)dS = 0 (5.1) 

Because of the symmetry assumption, the central inertia dyadics for F, S, and the 
overall system B can be written respectively as 

I f = If (b x b x + byb y ) + Jfb z b z (5.2) 

1$ = I* (b x b x + byby) + J s b z b z (5.3) 



and 
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I — I (b x b x + byby) + Jbjbj 


(5.4) 


We also have that 


I =I« + If + M« b 2 + nif a ^ 


(5.5) 


and 


J = J, + J f 


(5.6) 


Note also that the distances a and b of Fig. 1 are given respectively by 


a = 


Ma — c 

(M g + m f ) 


(5.7) 


b = 


mf 


(M, + m f ) 


(5.8) 


As propellant mass is lost, the system mass center B* moves forward toward S*. 
Combining (5), (7) and (8), we have 


I=I s + I f + pc 2 


(5.9) 


where 



■ I _ m f M s 

K (M s + m f ) 


(5.10) 


Letting the inertial angular velocity of S be 


co = 0) x b x + o>y by + ©I b z 


(5.11) 


the inertial angular acceleration is 


a = d) x b x + ©y by + d>z b z 


(5.12) 


(4) and (12) then give 


I • a = I (cb x b x + cby by) + J <b* b z 


(5.13) 


and (4), and (11) lead to 


to x I • CD = (J - 1) CDz (o) y b x - to* by) 


(5.14) 
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We also have from (4) and (11) that 



0) = i (oo* b x + tOy by) + j Gh b z 


(5.15) 


Finally, because of the simplifying assumption made so far. 


1 . 


p r x (g) x r) (v r n ) d5 = - m f [(4 + r iA) (g>x b x + ©y by) + (R 1 / 2 ) o>z bj (5. 16) 


where l t is the axial distance between F* and the nozzle exit plane. By combining (1) and 
(13)-(16), we have that the scalar equations of attitude motion for the system are 


I to* + ( J - 1 ) coy coz + [ i - m f ( l\ + R 1 / 4 )] a>x = 0 (5.17) 


I d)y - ( J -I ) to* Wz + [ i - riif ( 4 + R 1/4 )] COy = 0 (5.18) 


J cbz + ( Jf — rhf ) a>z = 0 (5.19) 


where I and J are given by (9) and (6) respectively, and [see (9)] 



80 


i = if + pc 2 + 2pcc 


( 5 . 20 ) 


j = jf 


( 5 . 21 ) 


5.2 Non-Dimensional Equations of Motion 


In this section, we render the equations of attitude motion very compact through the 
introduction of a set of dimensionless parameters. This allows us to perform further 
studies of these equations in a way that provides great insight into the most important 
factors that govern the behavior of the system. We start by introducing the non- 
dimensional time. For the system of Fig. 1, the instantaneous mass-of the cylindrical 
propellant is 


mf = mfo - m t 


( 5 . 22 ) 


where 


m =- m 


( 5 . 23 ) 



and is taken to be a positive constant that represents the rate of fuel consumption. We 
assume that t = 0 at ignition, when mf = rrifo, and that at burnout, t = t b and mf = 0. 
Hence, (22) gives, at burnout. 


tb = m ro/m 


We define a dimensionless time, x, as 


= X = jm_ t = at 
tb m fo 


Thus, x= 0 at ignition, and x = 1 at burnout. Also, the quantity cr 
time and is an appropriate scaling factor for time. 

We have from (25) that 



A. 

dx 


and the equations of motion, (17)-(19), can be re-written as 


(5.24) 


(5.25) 


has the dimension of 


(5.26) 
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. _ (j - I) COx COg - a [i - m' (le + R1/4)] 0)y 


COy = 


la 


(5.28) 


«z = 


- [j' - m' Rf/ 2 ] COg 

J 


(5.29) 


where the prime indicates derivative with respect to x , and the inertia properties I, J and 
mf are assumed to be expressed in terms of x . 

In the next phase of the non-dimensionalization process, we use mfo and R as scale 
factors for mass and length respectively, so that the scaling factor for inertia scalars 
becomes mm R 2 It is also convenient to introduce, at this point, the following non- 
dimensional parameters : 


Yi = / R , Y 2 = k sz / R , V|/ = M s / mo , 8 = L /R' 

P = Ri / R , 81 = ^1 / R , 82 = ^2 / R j 


(5.30) 


where k s and k sz are, respectively the transverse and axial radius of gyration of S. 

We now divide the numerator and denominator of the right hand side of (27) by 


mm R 2 to obtain 
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— 1 4 — coy CDz-a — 

mroR 2 mroR 2 / [mm R 2 mm R 2 

I ~ 


(l\ + R lAll to* 


If both sides of (31) are divided by a, we get 


Similarly, 




i-jl (4) 2 +2_ 


S y = (J/I-l)S x ^-l — 


J'_ JH- !— 
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where a>x, (Dy, C0z are non-dimensionalized angular velocity components, and I 
normalized inertia scalars. 


Following a strategy similar to that utilized in chapter 3, we let 


’ mpo 


fell 

\RI 



and 


6( T ) = j'_JL t 

w nifo 2 


Next, we multiply (32) by to* and (33) by CDy and add : 




(34) becomes 


■^-ah = -^-(dz 
dx j 


, J are 


(5.35) 


(5.36) 


(5.37) 


(5.38) 
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(37) and (38) lead to 


o) t = ©to exp 



\ 


i 



(5.39) 


and 


©z = ©zo exp - 


i 



(5.40) 


where 

©t = ^ (5.41) 


As in the case of the cylinder, the sign of the functions <|> and 6 determine whether 
the angular rates grow or decay with time. A positive sign implies decay and a negative 
sign signals divergence. 

From Fig. 1 


m f0 = pro it R 2 L 


(5.42) 
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where pro is the mass density of F. Continuity equation gives 


_ f R ‘ 

m = rhf = -m = 27t I p u(r) r dr 

Jo 


( 5 . 43 ) 


We recall here that p is the density of the fluid products of combustion at the exit plane. 
For a uniform velocity profile at the exit plane, u(t) is constant, so that 


m = riif = — m = — pjiRjU 


( 5 . 44 ) 


We then have an alternate expression for a as 


a .A. P" R i“ = _pJRi|>) = 11 b 2 i 

m ° pro ft R 2 L Pm'R' W ,P L 


( 5 . 45 ) 


where rj is the ratio p / pro. Note also that 



= a 5 


T1 


R 


( 5 . 46 ) 
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To make further progress with our analysis, we will now focus attention on two specific 
bum scenarios for F that are closest to what occurs in real rockets. 


5.3 Centrifugal Burn 

The fuel bums radially outwards in a symmetric manner (see Fig. 2) so that F* 
remains fixed in S; however, B* still moves toward S*. From Fig. 2, 

m f = p f 7tL(R 2 -r2) (5.47) 


U 



Fig. 5.2 


Centrifugal Bum 


From (47) and (44) 


^-Pf* L *2-p.« R ?» 


(5.48) 


Hence 




(5.49) 


and 


2 „ ^l 2 u 

‘ 2 =T 1 T 


t = T| 


R 2 u 
L a 


(5.50) 


From (50) and (45) 


r 2 P u a 

i— =n n x =“-T = x 

r 2 La a 


(5.51) 



and so. 


If 

m, 


from which we have 


Similarly, 


Jf 


and 


m f _ p f k L (R 2 - r 2 ) = l r 2 _ 1 ^ 
mro p f 7tLR 2 R 2 



4 12 12 4 


( 5 . 53 ) 



s 2 

12 2 


( 5 . 54 ) 


Jf _ mf (R 2 + r 2 ) _ y !+t - li Zl3. 
m„R 2 m ° 2R 2 ' 2 


( 5 . 55 ) 



- x 


( 5 . 56 ) 
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Because (see Fig. 2) 


we have 


and 


From (10) and (57) 


(J.C 2 

nv, R 2 


so that 


c = ^ 2 + L/2 


c = 0 


2 n c c’ = 0 


(m f /m 0 )(M s /m 0 ) I t 2 + Ll2 f_ (g ,£\2 

mf/mo + Ms/mo \ R I (l-x + y)* 2 2 / 


yL—tst+Ji.) 2 

m 0 R 2 (l-x + y) 2 ' 2/ 


(5.57) 


(5.58) 


(5.59) 


(5.60) 


(5.61) 



(20) yields 


-• -■ uc 2 2 u cc' 
I =I f + -£ + -^ 

m c R 2 m 0 R 2 


(5.62) 


We thus have, from (62), (54), (61) and (59) that 


r 11 

12 2 (\|f+l-x) 2 




) ; 


(5.63) 


Similarly, (21) and (56) give 



(5.64) 


Now (see Fig. 2), 


k 

R 




(5.65) 


and (52) leads to 
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mfo m ro 



(5.66) 


Hence, 


_dl 

mro 




(5.67) 


and (35), (63) and (67) yield 



V 2 

(\jf+ 1 -t) 2 




(5.68) 


In the same way, (36), (64) and (66) give 


0(t) = 



(5.69) 


5.3.1 Spin Rate Analysis 


From (69), the function 0(t) is a linear function of x ; it has a slope of -1 and is 
positive and equal to P 2 12 at x = 0, and is equal to P 2 1 2 - 1 at x = 1. Fig. 3 shows 0 
for various values of (3 . 
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Fig. 5.3 The function 0 for centrifugal bum 


We conclude from the figure that the spin rate always decreases at the beginning of the 
bum. If (3 > fl, the spin rate continues to decrease all the way to burnout For p < fl, the 
spin rate bottoms out during the bum, and then increases during the later phases of the 
bum. The turning point for this latter case occurs at 


e=- x + P 2 h = Q 


(5.69) 


That is, at 


x = 



(5.70) 



This result is consistent with, and augments that obtained for the cylinder, for which p = 1 
[ see (4.57) and Fig. 4.8]. The same result is validated by Fig. 4, obtained by numerical 
integration of (40). 


3 



Fig. 5.4 Centrifugal bum : the effect of expansion ratio p on spin rate 


5.3.2 Transverse Rate Analysis 

From (68), the slope of the function 0(t) is 


0 = — *- 
v 2 


2\|T 


(\|/+1-t) : 


■(«* 




(5.71) 


and remains negative throughout the bum. 
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At ignition ( x = 0 ), 



\|/ 2 

(y+1) 2 




(5.72) 


and at burnout ( x = 1 ), 


<t>i = 




- , P 

+ 5i 5 + 

4 



(5.73) 


It is clear from (72) and (73) that (j>o > 4>i» and that it is possible for any of these quantities 

f 

to be negative or positive. These facts, coupled with the fact that $ is negative, mean that 
the three scenarios depicted in Fig. 5 are possible, depending on the values assigned to the 
parameters. 



Fig. 5.5 The function <J> for centrifugal bum 
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We conclude from Fig. 5 that the magnitude of the transverse angular velocity for 
centrifugal bum can (i) decrease from ignition to burnout, (ii) decrease initially, then 
diverge later, or (iii) increase from ignition to burnout, depending on system parameters. 
Factors such as the distance of the exit plane from the combustion chamber (5i), shape 
factor of the propellant grain ( S), location of the propellant grain ( 82 ), nozzle expansion 
ratio (P), and the relative amount of propellant in the system (\p), all affect the transverse 
angular velocity. The main difference between the general axisymmetric system being 
considered here and the variable mass cylinder considered in chapter 4 is that case (iii) in 
Fig. 5 does not exist for cylinders. In other words, the transverse angular velocity of a 
variable mass cylinder always decreases initially, while a general axisymmetric system can 
exhibit divergent transverse rate throughout the propellant bum. Fig. 6 shows the 
transverse angular velocity magnitude as obtained by numerical integration of (39). 


X 



Fig. 5.6 Centrifugal bum : the effect of various 5 on cross-spin 
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5.4 Uniform Burn 


We recall that in uniform bum, the dimensions of F remain unchanged, while its 
density decreases continuously with time - see Fig. 7. In this case, mf has the same 
expression as mro in (42). This expression and (44) lead to 




^=p f nLR 2 = -p7tR 1 2 u 


(5.74) 


Hence 


‘Pf 

d(p f ) = -p 
pro 



(5.75) 
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so 


Pf = PfO~ P 


R^u_ t 

R 2 L 


From (74) and (76) 


mf = pro n R 2 L - p n R* u t 


mro 


Pro n R 2 L - p k Rj 2 u t 
ProJiLR 2 


= 1 -at = 1 -x 


and so. 


t, — ^ — a 


m„ R 2 


m n 


_ELl_ + _Li_ 


14 R 2 12 R 2 


=(>-441 


from which we have 



4 12 


J Jf _ m f R 2 _ 1 zl 
m 0 R 2 m ° 2 R 2 2 


(5.76) 


(5.77) 


(5.78) 


(5.79) 


(5.80) 


(5.81) 
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and 



(5.82) 


As before (see Fig. 7), 


and 


From (10) and (83) 


c = ^ 2 + L /2 


c =0 


2 \i c c' = 0 


(5.83) 


(5.84) 


(5.85) 


lie 2 _ (mf / m 0 ) (Mg / m 0 ) | ^ 2 + L/2| 2 _ (1~*)¥ [g . s\ 2 
m c R 2 mf / mo + M s / m 0 V R / (l-x + y)* 2 2* 


(5.86) 


and 


m 0 R 2 (l-x + v) 2 ' 2/ 


(5.87) 


We thus have, from (62), (80), (85) and (87) that 
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r=_&i_i_ 11 

2 (y+i_t) 2 


M 


(5.88) 


Similarly, (21) and (82) give 


J ’ = - J f = - 1 12 


(5.89) 


Now (see Fig. 7), 


and (78) leads to 


4 = £l + _l_ = Si+ & 

R R 7 R ^7 


Hence, 


m' _ m f 
mfo mfo 



(5.90) 


(5.91) 


jnl 

mro 




(5.92) 


and (35), (88) and (92) yield 



¥ 2 

(\|/+ 1 -t) 2 




(5.93) 


In the same way, (36), (89) and (91) give 



5.4.1 


Spin Rate Analysis 


From (94), the function 0(t) is a constant. Fig. 8 shows 0 for various values of 



Fig. 5.8 The function 0 for uniform bum 


If P > 1, the spin continues to decrease all the way to burnout. For P = 1, the spin rate 
always equals its initial value. For P < 1, the spin rate increases from the beginning to the 
end of the bum. The same result is validated by Fig. 9, obtained by numerical integration 
of (40) 



102 


S 

(Vi 





>K = 2 

a = 0.01 
5 s 10 
5 t = 2 
&2 = 3 
Yt * 1-2 


Fig. 5.9 Uniform Bum : the effect of expansion ratio (3 on spin rate 

5.4.2 Transverse Rate Analysis 

From (93), the slope of the function 4(x) is 


♦ - + &rY 

(w+l-x) 3 ' 2* 


(\|T + 1 — l) 


and stays negative throughout the bum. At ignition ( x = 0 ), 


^ = K + V + 6l5 + V)"4-[^i( 82+ 2 ) 


(Y+1V 


(5.95) 


(5.96) 




and at burnout ( x = 1 ), 



(5.97) 


Hence <{>o > <|>i, always, and either of these quantities can be negative or positive. Since we 

f 

also have that $ is always negative, we have the same scenario obtained for centrifugal 
bum. Numerical simulation results provided in Fig. 10 below confirm these inferences. 


a 

lO 

cm 



y = 15 

P-1 

a = 0.01 

81 = 3 

82 = 2 
Yi = 1-2 
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CHAPTER 6 

FORCED MOTION 
OF AXISYMMETRIC SYSTEMS 


All the development of the last two chapters have been restricted to systems 
subjected to zero external torque. The only forcing function present has been a perfectly 
axial thrust, whose line of action passes through the system mass center at all times, and 
thus produces no couple on the system. Great pains are taken in the design of real rockets 
to minimize thrust misalignments. Yet, it is unrealistic to expect that die system mass 
center, which generally moves during rocket bum, will remain positioned on the line of 
action of the resultant thrust vector throughout the propellant bum. Thus, torques due to 
thrust misalignment are often inevitable. There are also other potential sources of 
extraneous torques on a rocket system in flight, including aerodynamic and gravitational 
moments. For this reason, we will, in this chapter, develop relationships that describe the 
effects of externally applied moments on the attitude behavior of the same axisymmetric 
variable mass system that was studied in the previous chapter. 
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6.1 Equations of Attitude Motion 

The equations of forced attitude motion of the system of Fig. 5.1 are really the same 
as (5.17)-(5. 19), except that the right hand side of each of these equation is non-zero. 

Thus, we have 


Ici)x + ( J - 1 )©y ©z + [ i-m f ( ti + R 1/4 ) ] ©X = M x (6.1) 


I ci)y -( J-I ) ©x ©z + [ I - ihf ( tl + R1/4 )] ©y = My (6.2) 


J©z + ( jf-m f R?/2)©z = M z (6.3) 


The forcing functions M* , My and M* are, in general, functions of time. This is so 
because at least the contributions from thrust misalignments vary with time since the system 
mass center moves relative to the system base structure. 

We start the solution process by letting 


\ 2 = (if— rhfRi 2 / 2) / J 


(6.4) 


and 


M z = M z / j 


(6.5) 



(3) then becomes 


^- + Mt)<Oz = M z 


( 6 . 6 ) 


The general solution of (6) is 



c _ 1 r \ 

/ 

CDz = 

Oho + I Mjj(s)exp I X^x) dx ds 

exp - I A. 2 (t) dr 


Jo \ Jo ) 

\ Jo / 


(6.7) 


This solution can be thought of as the sum of two components : a first term that represents 
die homogenous part due to initial condition, and a second term - the particular solution - 
that is contributed by the forcing function. If J, j f , riif and M z are known functions of 
time, and the initial spin rate is known, (7) determines uniquely the spin rate (Oz of the body 
S. Note that Ot)z depends only on X 2 and M z , and is independent of M x and M y . This is 
true because of our assumption of axisymmetry. 

We are now ready to investigate the transverse angular velocity components of the 
system. To this end, we introduce the following complex quantities : 


+ j ®y 


( 6 . 8 ) 


M X y — M x + j My 


(6.9) 
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where 


j = V^T 


We also define 


Xi = [i - mi (4 + R 1 / 4 )] /i 


o = ( J-I)/i 


and 


M X y = M xy / J 


Adding j times (2) to (1), we obtain 


I QJxy - j (J — I) Q>z (j COy+CD^ + Xl <0 X y = M X y 


or 


( 6 . 10 ) 


( 6 . 11 ) 


( 6 . 12 ) 


(6.13) 


(6.14) 
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^ + P 2 (t)o) xy = M xy 


where 


P 2 <t) = Xi — a coz 


The solution to (16) is 



r _ If* ) 


If 1 

CDxy - 

CDxyo+l M xy (s) exp I p 2 (t) dx 

ds 

exp - I P 2 (x)dx 


Jo [Jo ) 


\ Jo , 


or equivalently, 


(Ox +j (Dy = 


| 0^0 + j CDyO + 


[ 


[M X (s)/l+j My(s)/ l] 


(exp[r (s)]} [cos0(s)- j sin©(s)] ds} (exp [- T (s)]} [cos0(t) + j sin0(t)] 


where 


r(t)= 



Xi dr 


(6.15) 


(6.16) 


(6.17) 


(6.18) 


(6.19) 


and 



0(t)= 



(a ci)z) dx 


( 6 . 20 ) 


The functions r(t) and 0(t), which have fundamentally different properties, govern the 
character of the transverse motion of the variable mass system under consideration. r(t) 
determines the amplitude of the angular motion and 0(t) affects motion frequency. Explicit 
formula for to* and Ct>y can be obtained from (1 8) by separating the real and imaginary parts 
of the expression; they are : 


<Dx(t) = e - T ( l ) [cDxo cos0(t) - Q)yo sin6(t)] + 

[ e ^ r(s | - ~ r -^ (m, (s) cos[e(s) - e(tjj ♦ My (s) sm[e(s) - e(t)]) ds 

/to 


( 6 . 21 ) 


and 


(Oy(t) = e - ^ [tt>xo sin0(t) + (Dyo cos0(t)] + 

f EM j My (s) cos[0(s) - 0(t)] - M x (s) sin[@(s) - 0(tj] ) ds 

/to 


( 6 . 22 ) 


We observe that in general the spin torque, M z has no effect on the transverse angular 
velocity, nor is the spin rate affected by any of the transverse torques. 
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6.2 Change of the Independent Variable 

We now introduce a new independent variable in order to simplify somewhat the 
equations of attitude motion and their solution. We exploit here, a transformation defined 
earlier and first proposed by Jarmolow (1957) : 

0(0= ( oh^f^-d £=( CDz ad£ (6.23) 

Jo Jo 


0, has the dimension of angle (radians) and is related to the spin angle. It will be used as 
the new measure of time. 

Note that 


d@ = G)z 


I 


(6.24) 


SO 



(J*-I*) d 

I* d© 


(6.25) 


The starred quantities are assumed to be expressed in terms of the new variable 0. This is 
indeed possible because (23) gives 0 as a function of t, and can be inverted to yield the 



expression of t as a function of 6. Equations (1) and (2) can be re-written in terms of @ as 
follows : 
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(J* - 1* ) ^ + ( J* - 1* ) a* <* + k - ! *2 l (/; 2 + Rf/J ^ = K (6.26) 


d0 d0 


a^(r _ J* - 1 * ) ©x ©z + [^-^(fe 2 + Rl/J ©z ©y = M*y ( 6 . 27 ) 

d0 Ld0 d0 Jr 


Dividing through by (j* - 1*) 0 ) 2 , one obtains 


d0 Ld© d0 Jl ©z(j*-I*) 


(6.28) 


sat-^+fdt-asHtf t Riu|aL ( &29> 


i* ©z(r-r) 


Then, multiplying (29) by j and adding to (28), we have 


d© X y (, * . 


(Xi — j) ©xy — Mxy 


( 6 . 30 ) 



where 


to xy (e) = Q) x ( 0 ) + jo)y( 0 ) 


(6.31) 


X 


* 

1 




(6.32) 


and 


M xy = 


j My 

COz (j* - I*) 


(6.33) 


If we define T ( 0 ) as 



(6.34) 


then the solution of (30) for co* and coy in terms of the new variable 0 can be written 
explicitly as 



(Ox ( 0 ) = (cOxo cos0 - (Oyo sin0) e _ ^ (®) + 

— — [m^ (s) cos(s - 0 ) + (s) sin(s - ©)] ds 

(0zU*-n 


(6.35) 
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Cl)y (0) = [tDxO sin© + COyO cos©] e~ I' ’(©)+ 


r 

JO 


— — [m£ (s) cos(s - ©) - MJ (s) sin(s - ©)] ds 
tOz U - 1 ) 


(6.36) 


In order to obtain explicit expressions for the known functions such as I, J, M x and 
My, the function t(0) must be determined. This may or may not be a trivial task. The 
integrand of (23) is a rational function of time whose integral may not have a closed form 
analytical expression. However, a numerical solution is always possible, so we can be 
sure that the complete solution sketched above can always be obtained by numerical means 
- at least 


6.3 Special Case - Rigid Body Motion 

The above solution can be used to investigate the attitude motion of an axisymmetric 
rigid body simply by setting 


sUl = sUl -dmf-p 
d© d© d© 

(32) and (34) then yield 


(6.37) 


r*(©) = 0 (6.38) 


and (35) and (36) simplify as follows 
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(Ox ( 0 ) = (cOxo cos© - <»yo sin©) + 

J - -J ^ [m*x(s) cos(s - @) + (s) sin(s - ©)] ds 

0 (OzU -I ) 


(6.39) 


<Oy (©) = [g>xo sin© + (Dyo cos©] 


(6.40) 


re 


J ! -v ,. 1 cos(s - ©)- (s) sin(s - ©)] ds 

0 COz(J -I ) 


In this form, the physical significance of the various parameters of the problem is not so 
obvious. Simpler and more revealing expressions for (Ox and G)y are obtained for the case 
of constant body-fixed torques. For instance, if we take M x and M y to be constants, and 
M z = 0, (39) and (40) further simplify to 


o>x(0) = [(Ox oCO s(©)- (Oyo sin(©)] + 

(j* - 1*) (Oz -M ^ 1_ cos ( 0 SI 


(6.41) 


(Oy(0) = [(Oyo cos(0) + (Oxo sin!©)] + 
^ {m£ sin(0) + [l- CO 


(6.42) 


where 
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© (0 = —j— ■ COzO t 


(6.43) 


We note here that the spin rate remains constant at its initial value, as can be seen 
from (4), (5) and (7). (41) and (42) indicate that transverse velocities increase with 
transverse torques and decrease with the inertia difference, J* - 1*, and initial spin rate. 

Spinning a rigid body about one of its principal axes provides the so-called spin 
rigidity or gyroscopic stiffness; that is, a resistance to external disturbances during motion. 
This is the well known advantage of spin stabilization. In a spin stabilized system, the 
transverse angular rates to* and CDy are usually the result of unforeseen disturbance, and 
are generally undesirable. If M x and M y are viewed as constants, body-fixed disturbance 
torques in (41) and (42) above, then, these equations show clearly that high spin rates will 
reduce the effects of the disturbance torques on vehicle motion. This indicates that the 
nominal spin rate of a spin stabilized vehicle should be given the highest value compatible 
with the vehicle’s structural integrity. A fact that agrees with physical intuition. It can thus 
be said that as the spin rate increases, the stability of a given system increases as measured 
by the reduction in angular rates introduced by a given disturbance. 

For a single spinning body, (41) and (42) also show that gyroscopic stiffness 
depends in a crucial way on the system inertia property, J* - 1*, or the difference between 
the spin and the transverse moment of inertia of the body. Spin rigidity vanishes as the 
inertia difference tends to zero. In other words, bodies with spherical symmetry can have 
no gyroscopic stiffness, irrespective of the value the spin rate. On the other hand, heavily 
prolate or oblate bodies are seen to be quite stable in spin. 



CHAPTER 7 


CONCLUSION 


The study presented in this report deals with the dynamic behavior of spinning 
bodies that lose mass while in modon. The study began with the determination of the 
complete equations of such systems using one of the most powerful methods of analytical 
dynamics — Kane’s formalism. The resulting equations arc then applied to the study of 
the attitude dynamics of space-based variable mass systems using mathematical models that 
vary from a simple cylinder to general axisymmetric systems. The emphasis of these 
efforts was on information extraction. Special mathematical techniques were used to 
extract as much information as possible about the motion of the system, without actually 
attempting a full solution of the equations of motion. 

In the first example presented, a variable mass rocket system is represented by a 
simple cylinder with four possible mass loss scenarios. Results obtained indicate that such 
a system can fly stably in a torque free environment if the ratio of the radius to the length is 
small. Large and short cylindrical vehicles can have stability problems, especially if they 
are subjected to radial bum. This study goes far beyond previous work in this area by 
clearly identifying the exact inertia conditions and the precise stages of vehicle motion at 



which the character of motion begins to change. For cases where there is stability issue, 
limiting values of motion parameters are determined 

Attention was next turned to a slightly more complicated model — a general 
axisymmetric variable mass system. Although important differences were found in the 
details of the motion of the general axisymmetric system as compared to that of the simple 
cylinder, the main stability results did not change. It was found in this second example, 
that in addition to the shape factor of the system propellant, the location of the propellant 
grain within the system as well as the nozzle expansion ratio all have significant effects on 
the attitude motion of a rocket. 

In the last section of this study, the combined effect of mass loss and external 
transverse torques on system motion is examined and solutions to the equations of motion 
in this case are presented Overall, this study augments considerably, previous knowledge 
of the dynamics of variable mass systems and gives great insight into the effects of mass 
variation on the attitude behavior of spinning bodies. 

This study shows precisely how the system’s attitude response is tied to such 
system parameters as nozzle geometry, combustion chamber geometry, location and 
relative mass of the propellant In this sense, the study will be of great use to designers in 
that it indicates what must be done to avoid motion instability of the type witnessed in the 
80’s on Star 48 propelled space vehicles. 
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